Non Gaussian extrema counts for CMB maps 
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In the context of the geometrical analysis of weakly non Gaussian CMB maps, the 2D differential 
extrema counts as functions of the excursion set threshold is derived from the full moments expansion 
of the joint probability distribution of an isotropic random field, its gradient and invariants of the 
Hessian. Analytic expressions for these counts are given to second order in the non Gaussian 
correction, while a Monte Carlo method to compute them to arbitrary order is presented. Matching 
count statistics to these estimators is illustrated on fiducial non Gaussian "Planck" data. 
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Random fields are ubiquitous phenomena in physics 
appearing in areas ranging from turbulence to the land- 
scape of string theories. In cosmology, the sky-maps of 
the polarized Cosmic Microwave Background (CMB) ra- 
diation - a focal topic of current research ~ is a prime 
example of such 2D random fields. Modern view of the 
cosmos, developed primarily through statistical analysis 
of these fields, points to a Universe that is statistically 
homogeneous and isotropic with a hierarchy of structures 
arising from small Gaussian fiuctuations of quantum ori- 
gin. While the Gaussian limit provides the fundamental 
starting point in the study of random fields [IHS], non- 
Gaussian features of the CMB fields are of great interest. 
Indeed, CMB inherits a high level of gaussianity from ini- 
tial fluctuations, but small non-Gaussian deviations may 
provide a unique window into the details of processes in 
the early Universe. The search for the best methods to 
analyze non-Gaussian random fields is ongoing. 

In paper [4] the general invariant based formalism for 
computing topological and geometrical characteristics of 
non Gaussian fields was presented. The general formulae 
for the Euler characteristics to all orders has been de- 
rived, which encompasses the well known first correction 
[5] and which was later confirmed to the next order by [B] . 
We now focus on the statistics of the density of extremal 
points which follows directly from the formalism of [1]. 
The goal of this paper is to provide an explicit recipe 
on how to use this formalism in practice on idealised 2D 
CMB "Planck" -hkc data. 



I. EXTREMA COUNTS 

Extrema counts, especially that of the maxima of the 
field, have long application to cosmology [e.g. [3], however 
theoretical development have been mostly restricted to 
the Gaussian fields. The statistics of extrema counts, as 



well as of the Euler number, requires the knowledge of 
the one-point JPDF P{x,Xi,Xij) of the field x, its first, 
Xi, and second, Xij, derivatives |llj . Extrema density is 
an intrinsically isotropic statistics given by [U [7] 
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Under the condition of statistical isotropy of the field, the 
essential form for the JPDF is therefore given in terms 
of the rotation invariants — x itself, the square of the 



magnitude of the gradient q 



X2 

^ of the Hessian 



and the two 



invariants Ji = Ai -)~ A2, J2 = (Ai — A2 
matrix xij (where A^ are the eigenvalues of the Hessian) . 
Introducing ( = (x + "fJi)/\/l — 7^ (where the spectral 
parameter 7 = —{xJi) characterizes the shape of the 
underlying power spectrum), leads to the following JPDF 
for the Gaussian 2D field 



tr2D = :z- cxp 

ZTT 






The invariant form for the extrema counts 



driest 
dx 



dJidJ2 



8^2 yr 



: exp 



■T 



-A - J 2 



(2) 



\Ji-J2 



then readily recovers the classical results [H [31 [7] when 
the limits of integration that define the extrema type 
are implemented, namely Ji G [— cx),0], J2 G [0, Jf ] for 
maxima, Ji € [0, 00], J2 G [Oj «^i ] for rninima and Ji e 
[— oo,cx)],J2 G [Jf,(X)] for saddle points. 

In [4] we have observed that for non-Gaussian JPDF 
the invariant approach immediately suggests a Gram- 
Charlier expansion in terms of the orthogonal polyno- 
mials defined by the kernel G2D- Since C,, g^, Ji and 
J2 are uncorrelated variables in the Gaussian limit, the 
resulting expansion is 
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where terms are sorted in the order of the field power n 
and '^l ■ j! ; Jg ~" stands for summation over all combi- 
nations of non-negative i, j, fc, I such that i + 2j + k + 21 
adds to the order of the expansion term n. 

The Gram-Charlier coefficients, ( C q^^ ■h'" J2) = 

\ /GC 

i-iy+'jm (H, (C) L, (g2) Hk ( Ji) U ( J2))^ that appear 
in the expansion can be related to the more famil- 
iar cumulants of the field and its derivatives (we use 
( )ni for statistical moments while reserving ( ) for 
statistical cumulants), actually being identical to them 
for the first three orders n = 3,4,5. Lookup ta- 
bles of the relationship between Gram-Charlier cu- 
mulants and statistical cumulants can be found at 
'http : //www . lap . f r/users/pichon/Grgmi/I As an illus- 
tration, one sixth order non trivial cumulant would be 

( Jf J2C) CG = ( Jf J2C) + (J!) {J2O + 3 ( Ji J2) ( J?C). It 
is prudent to stress that the Gram-Charlier series expan- 
sion is distinct from the perturbative expansions. For 
instance, while the linear Edgeworth or /nl expansion 
match solely to the first order n = 3 Gram-Charlier coef- 
ficients, quadratic terms require knowledge of the Gram- 
Charlier terms to n = 6, while the cubic ones to n = 9. 
Integrals over Ji and J2 for extremal points can be car- 



ried out analytically even for the general expression pi. 
Different types of critical points can be evaluated sepa- 
rately by restraining the integration domain in the J1-J2 
plane to ensure the appropriate signs for the eigenvalues. 

The effect of the non-Gaussian cubic correction on the 
total number of the extrema of different types is given by 



'^max/niin 
^sad 



± 



1 

8V37ri?,2 

1 
4V37ri?*2 ' 



18((?2ji)-5(jf)+6(JiJ2) 



MttV^^R* 



(4) 



where we have restored (see note [11]) the dimensional 
scaling with R^ = cr 1/(72 , the characteristic separation 
scale between extrema. The total number of saddles, as 
well as of all the extremal points, rimax + JT-min + ^T-sad, 
are preserved in the first order (the latter following for 
the former, as topological considerations imply n^ax ~ 
"sad + "mill = coust), but the Symmetry between the 
minima and the maxima is broken. 

The differential number counts with respect to the ex- 
cursion threshold v are given by 
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where i^i, K^., K^ are polynomials with coefficients expressed in terms of the cumulants. Here we give explicit expres- 
sions for the first non-Gaussian order, while the next order can be found at the above mentioned URL. 

The term Ki{v^ 7) has a special role determining the Eulcr number x{^) via dx/ dv = d jdv (nmax + J^min ~ "-sad) ~ 
^2/7rexp(— 1/2/2)^1 (i^, 7). As such, its full expansion has been given in [4 , Eq. (7), and confirmed to the second 
order in [6]. To the leading non-Gaussian order 
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Introducing scaled Hermite polynomials ^-[^{v^a) = cr^"il„ (i'/ct), the polynomial K2{i^,j), the only one that 
determines the distribution of saddle points, can be written as 
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The remaining term, i^3(z^, 7) is the most comphcated one. It is expressed as the expansion in '}i^{v, \J\ — 7^): 
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Eqs (l5])-(l6| (together with the next order expansion avail- 
able online) are the main theoretical result of this paper. 



II. IMPLEMENTATION 

Evaluating these estimators requires computing the cu- 
mulants appearing in Eqs. ([7|)-([9]). In non-Gaussian mod- 
els where the field is represented by the functional of 
a Gaussian field this may be possible directly, while in 
general, as shown in [B], such cumulants can be found 
as weighted marginals of the underlying bispectrum, (to 
third order), trispectrum (to fourth order), etc.. On a 
sphere, the high order marginals are particularly cum- 
bersome and time consuming to compute, as they also 
involve the contractions of n — j Wigner symbols. Here 
we suggest a different route, based on the assumption 
that scientists interested in fitting extrema counts to non- 
Gaussian maps are typically in a position to generate 
realizations of such maps. In that case, it becomes rela- 
tively straightforward to draw samples of such maps, and 
estimate the corresponding cumulants. The HEALPix [5] 
library provides in fact a direct estimate of the deriva- 
tives of such maps up to second order, which is all that 
is required to compute the cumulants of the JPDF. 

As an illustration, let us generate sets of parameter- 
ized non-Gaussian maps using the package sky-ng-sim 
[S] of HEALPix. In this so called harmonic model, 
the PDF of the pixel temperature, T is given by 
exp(-r2/2cr^) II]^=o"i^«-^i(^/^o)|^ where d are nor- 
malization constants. In this paper, we use nside=2048, 
^max = 4096, n = 2, CTo = 1, ao = and vary a\ and a2- 
We also consider the second option of sky-ng-sim which 
produces non Gaussian field as even power, /? of unit vari- 
ance zero mean Gaussian fields. For each set of maps, we 
compute its derivatives, and arithmetically average the 
corresponding cumulants, using a code, map2cum relying 
on the HEALPix routine alm2map_der. Invariant vari- 
ables J\ and Ji on a sphere are defined via the mixed 
tensor of covariant derivatives J\ — x-^i'^ and J2 = KC;i'"' ■ 
The differential counts are then evaluated for a range of 
threshold, v S [—5, 5]. For each of these maps, the num- 



ber of extrema is computed by the procedure map2ext 
which implements the following algorithm: for every pixel 
a segment of quadratic surface is fit in the tangent plane 
based on the temperature values at the pixel of origin and 
its HEALPix neighbours. The position of the extremum 
of this quadratic, its height and its Hessian are computed. 
The extremum is counted into the tally of the type de- 
termined by its Hessian if its position falls within the 
original pixel. Several additional checks are performed 
to preclude registering extrema in the neighbouring pix- 
els and minimize missing extrema due to jumps in the fit 
parameters as region shifts to the next pixel. Masks are 
treated by not considering pixels next to the mask bound- 
ary. Pixel-pixel noise covariance can be included while 
doing the local fit. On noise- free maps the procedure 
performs with better than 1% accuracy when the map 
is smoothed with Gaussian filter with FWHM exceeding 
6 pixels. Both map2cuni and map2ext are available upon 
request. Figure [l] illustrates the very good agreement be- 
tween the theoretical expectation of the differential num- 
ber counts to the measured ones for both the harmonic 
and the power-law models. 

An alternative numerical procedure, which is likely to 
be more practical for expansion beyond the fourth order 
was also successfully explored for 2D topological invari- 
ants. Starting from Eq. p|, we re-express both the poly- 
nomials in Ji, J2, C, and g^ and G2D in terms of the six 
field variables, {x,Xi,Xij). We then construct formally 
the marginal G'^(x = (a^ii, 2^121 2;22)|a; = v,xi — X2 — 0), 
where the latter condition corresponds to imposing that 
we are seeking extrema of the field. It becomes straight- 
forward to draw large sets of 3 random numbers satisfy- 
ing Gy. For each triplets, x, and a given numerical set 
of cumulants, we then compute the argument, I(x) of 
the square bracket in Eq. (|3| (up to some given order), 
together with the two eigenvalues of the Hessian. For 
maxima (resp. minima, resp. saddle points), we replace 
I by if the two eigenvalues are not negative (resp. pos- 
itive, resp. of different sign). The sum over all triplets 
yields a Monte Carlo estimate of dnc-,^t/dv. The accuracy 
of the estimate depends on the extent of rejection while 
applying the extremal condition. 
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FIG. 1: Top panel: the predicted {solid line) number of maxima, saddle, and minima in Av — 0.25 bins as a function of 
the threshold, u, on top of the measured count from a single realization full-sky nside=2048 HEALPix map (histogram). 
The temperature field is smoothed with the Gaussian filter of 10 arcmin FWHM, resulting in R, « 5.5 arcmin « 3 pixels. 
The dashed line corresponds to the Gaussian prediction. The left panel corresponds to the Harmonic Oscillator model of non 
Gaussianity with ai — 0.6, a2 ~ 0.6, while the right panel corresponds to the power law non Gaussianity with j3 = 2. Bottom 
panel: the departure from Gaussianity for these two models as predicted (solid line) and measured (dashed line) for maxima 
(green), minima (blue) and saddle points (red). Note that the corrections of Eqs (|5|-(|6| (solid line) give a very accurate match 
to the measured PDF. As is seen, different models of non-Gaussianity can be distinguished by their effects on extrema. 



Note in closing that all the presented analysis is 
straightforwardly generalized to 3D (noticeably the 
Monte Carlo method), as shown in [TOj, to describe the 
large scale distribution of matter. Indeed in this context, 
the gravitational instability that nonlinearly maps the 
initial Gaussian inhomogeneities in matter density into 
the LSS, induces strong non-Gaussian features culminat- 
ing in the formation of collapsed, self-gravitating objects 
such as galaxies and clusters of galaxies. 
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